Project Report

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(plotly)

## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

library(shiny)
library(flexdashboard)
library(viridis)

## Loading required package: viridisLite

data_clean <- read.csv("data_final.csv")

Project Name: Examining Demographic Patterns in NYC Shooting Data

Project Members: Chenhui Yan (CY2772), Zhaokun Lin (ZL3544), Mingyin Wang (AW3693), Zebang Zhang (ZZ3309)

Topic

In this exploratory study, we aim to analyze the geographic and temporal distribution of shootings in New York City, with an emphasis on how socioeconomic factors influence these patterns.

Motivation

Understanding the dynamics of shooting incidents in New York City is crucial for enhancing public safety and fostering resilient communities. By examining demographic patterns in shooting crime data alongside datasets on high school graduation rates and poverty levels, our project aims to identify the socioeconomic factors that influence gun violence. Analyzing how disparities in education and income relate to shooting incidents provides valuable insights for targeted interventions.

Introduction

Gun violence is a significant public health and social issue New York City, where diverse communities and complex socioeconomic conditions create a unique landscape for understanding its contributing factors. This project examines the geographic, temporal, and socioeconomic patterns of shootings in NYC to identify critical influences on gun violence. We explore how shooting incidents are distributed across boroughs and neighborhoods and examine temporal trends, such as seasonal, monthly, and time-of-day variations. Additionally, we investigate the relationship between education and gun violence, focusing on whether lower high school graduation rates correlate with higher shooting prevalence. Finally, we analyze the connection between poverty and gun violence to determine if neighborhoods with higher poverty levels experience increased incidents of shootings.

Significance

This research will aid in informing targeted public health and safety interventions by identifying the communities most impacted by shootings. Understanding these patterns will enable policymakers to implement more effective violence prevention measures and allocate resources where they are needed the most.

Research Questions

Our project focuses on examining the socioeconomic and temporal factors that influence gun violence in New York City. Specifically, we aim to answer the following questions:

Geographic and Temporal Patterns of Shootings in NYC: - How are shooting incidents distributed geographically across different boroughs and neighborhoods in New York City? - What are the temporal patterns of shootings (e.g., month, season, time of day) in different areas of NYC?

Education and Gun Violence: - Is there an association between high school graduation rates and the prevalence of shootings across NYC neighborhoods? - Are neighborhoods with lower high school graduation rates more likely to experience higher rates of gun violence?

Poverty and Gun Violence: - How does the percentage of people living below the poverty line relate to shooting incidents across neighborhoods in NYC? - Are higher poverty levels associated with an increased frequency of shootings?

Data source

Original data

This report analyzes shooting incidents in NYC. You can find the original data before merging here. This data is from NYCOpenData, consists some basic information: occur time, coordinate of each incident, victim’s information, .etc.

data for additional information

• NTA shape data

Neighborhood Tabulation Areas (NTAs) are medium-sized statistical geographies for reporting Decennial Census and American Community Survey. NTAs were delineated with the need for both geographic specificity and statistical reliability in mind. Shapefile of NTA is get from NYC Planning website, you can find the data here, or you can download

• CDTA data

The Department of City Planning (DCP) created Community District Tabulation Areas (CDTAs) to closely approximate the 59 Community Districts of New York City for the purpose of reporting American Community Survey (ACS) data. You can find the data here, or you can download

• Borough Boundaries data

The Borough Boundaries dataset from the NYC Open Data portal provides detailed information about the geographical boundaries of the five boroughs of New York City: Bronx, Brooklyn, Manhattan, Queens, and Staten Island. You can find the data here

• Neighborhood data

Here is neighbourhood list for geo filter. Sourced from city or open source GIS files. Here is GeoJSON file of neighbourhoods of the city.

• Education data

The data source for education is shown here.

• Poverty data

The data source for poverty is here

Data processing

The original data only consists of some basic information and doesn’t show NTA, CDTA, Borough and neighborhood information of each incident, so we added these information to the original data base on the coordinate of each incident, we also add some demographic information: population, education level, poverty level, .etc to the data.

Here’s how we added the information based on all the data we have:

• enrich NYC shooting incident data with neighborhood-level information by performing spatial operations, converts the shooting data into a spatial format, joins it with neighborhood boundaries, and adds additional neighborhood attributes for deeper spatial analysis.
  detailed code

  df=read_csv("./data/NYPD_Shooting_Incident_Data__Historic__20241119.csv")
  
  # Ensure coordinate columns are numeric
  df = df %>%
    mutate(
      X_COORD_CD = as.numeric(X_COORD_CD),
      Y_COORD_CD = as.numeric(Y_COORD_CD)
    )
  
  # Convert to sf object with CRS EPSG:2263
  df_sf = st_as_sf(df, coords = c("X_COORD_CD", "Y_COORD_CD"), crs = 2263)
  
  # Download and read the GeoJSON file
  download.file(
    url = "https://data.insideairbnb.com/united-states/ny/new-york-city/2024-09-04/visualisations/neighbourhoods.geojson",
    destfile = "neighbourhoods.geojson",
    mode = "wb"
  )
  neighbourhoods = st_read("neighbourhoods.geojson")
  
  # Check and transform CRS of neighborhoods
  if (st_crs(neighbourhoods)$epsg != 2263) {
    neighbourhoods = st_transform(neighbourhoods, crs = 2263)
  }
  
  # Perform the spatial join
  df_joined = st_join(df_sf, neighbourhoods, left = TRUE)
  
  # Add neighborhood information to the data frame
  df$Neighborhood = df_joined$neighbourhood
  # Download and use the neighborhoods CSV file
  download.file(
    url = "https://data.insideairbnb.com/united-states/ny/new-york-city/2024-09-04/visualisations/neighbourhoods.csv",
    destfile = "neighbourhoods.csv",
    mode = "wb"
  )
  neighborhoods_csv = read_csv("neighbourhoods.csv")
  
  # Merge additional attributes if necessary
  df = df %>%
    left_join(neighborhoods_csv, by = c("Neighborhood" = "neighbourhood"))

• enrich NYC shooting incident data with neighborhood-level information by performing spatial operations, converts the shooting data into a spatial format, joins it with neighborhood boundaries (both 2020 and 2010 shapefiles), and adds additional neighborhood attributes for deeper spatial analysis.
  detailed code

  df$X_COORD_CD <- as.numeric(df$X_COORD_CD)
  df$Y_COORD_CD <- as.numeric(df$Y_COORD_CD)
  
  # Convert to sf object
  df_sf <- st_as_sf(df, coords = c("X_COORD_CD", "Y_COORD_CD"), crs = 2263)
  
  # Read neighborhood shapefile
  nta <- st_read("./data/nynta2020_24d/nynta2020.shp")
  
  # Transform CRS if necessary
  nta <- st_transform(nta, crs = st_crs(df_sf))
  
  # Perform spatial join
  df_joined <- st_join(df_sf, nta, left = TRUE)
  
  # Add neighborhood information to your original data frame
  df$NTA <- df_joined$NTAName  # Adjust based on actual column name
  
  ## add 2010 nta data
  df$X_COORD_CD <- as.numeric(df$X_COORD_CD)
  df$Y_COORD_CD <- as.numeric(df$Y_COORD_CD)
  
  # Convert to sf object
  df_sf <- st_as_sf(df, coords = c("X_COORD_CD", "Y_COORD_CD"), crs = 2263)
  
  # Read neighborhood shapefile
  neighborhoods_2010 <- st_read("./data/nynta2010_24d/nynta2010.shp")
  
  # Transform CRS if necessary
  neighborhoods_2010 <- st_transform(neighborhoods_2010, crs = st_crs(df_sf))
  
  # Perform spatial join
  df_joined <- st_join(df_sf, neighborhoods_2010, left = TRUE)
  # Add neighborhood information to your original data frame
  df$NTA_2010 <- df_joined$NTAName  # Adjust based on actual column name

• using the holidayNYSE() function to determine which dates are public holidays
  detailed code

  df$OCCUR_DATE <- as.Date(df$OCCUR_DATE, format = "%m/%d/%Y")
  
  # Define the range of years, handling NA values
  years <- seq(
    year(min(df$OCCUR_DATE, na.rm = TRUE)),
    year(max(df$OCCUR_DATE, na.rm = TRUE)),
    by = 1
  )
  us_holidays = holidayNYSE(years)
  us_holidays = as.Date(us_holidays)
  # Add a new column indicating whether the date is a holiday, year and month information of a date
  df = df %>%
    mutate(
      Is_Holiday = OCCUR_DATE %in% us_holidays,
      Year = year(OCCUR_DATE),
      Month = month(OCCUR_DATE)
    )

• calculates dawn and dusk times for each unique date in the dataset. Using NYC’s coordinates (latitude and longitude), the getSunlightTimes() function generates dawn and dusk times for each day
  detailed code

  df$OCCUR_DATE <- as.character(df$OCCUR_DATE)
  df$OCCUR_TIME <- as.character(df$OCCUR_TIME)
  
  # Combine OCCUR_DATE and OCCUR_TIME into a single datetime string
  datetime_str <- paste(df$OCCUR_DATE, df$OCCUR_TIME)
  
  # Parse datetime using lubridate
  df$OCCUR_DATETIME <- ymd_hms(datetime_str, tz = "America/New_York")
  # Extract date from OCCUR_DATETIME
  df$DATE <- as.Date(df$OCCUR_DATETIME)
  
  # Get unique dates
  unique_dates <- unique(df$DATE)
  latitude <- 40.7128
  longitude <- -74.0060
  # Calculate dawn and dusk times
  sun_times <- getSunlightTimes(
    date = unique_dates,
    lat = latitude,
    lon = longitude,
    keep = c("dawn", "dusk"),
    tz = "America/New_York"
  )
  
  # Merge with the original dataframe
  df <- df %>%
    left_join(sun_times[, c("date", "dawn", "dusk")], by = c("DATE" = "date"))
  
  # Determine if the sky is dark considering twilight
  df <- df %>%
    mutate(
      Sky_Is_Dark = OCCUR_DATETIME < dawn | OCCUR_DATETIME >= dusk
    )
  df <- df %>%
    select(-DATE, -dawn, -dusk)

• reads in population data from an Excel file, filters relevant rows, selects specific columns, and merges this population data with the original dataset.
  detailed code

  popu <- suppressWarnings(
    read_excel("./data/nyc_detailed-race-and-ethnicity-data_2020_core-geographies.xlsx", 
               sheet = 1,
               range = "A4:I2599",
               col_types = c("numeric", "text", "text", "text", "numeric", "text", "text", "numeric", "numeric")) |>
      filter(`Orig Order` >= 2334 & `Orig Order` <= 2595) |>
      select(GeoName, NTAType, Pop) |>
      rename(Total_population = Pop)
  )
  
  df <- df |>
    left_join(popu, by = c("NTA" = "GeoName")) |>
    mutate(
      NTAType = case_when(
        NTAType == 0 ~ "Residential",
        NTAType == 9 ~ "Park",
        NTAType == 8 ~ "Airport",
        NTAType == 7 ~ "Cemetery",
        NTAType == 6 ~ "Other Special Areas",
        NTAType == 5 ~ "Rikers Island",
        is.na(NTAType) ~ "Unknown")
    )

• add CDTA data similar to NTA
  detailed code

  df$X_COORD_CD <- as.numeric(df$X_COORD_CD)
  df$Y_COORD_CD <- as.numeric(df$Y_COORD_CD)
  
  df_sf <- st_as_sf(df, coords = c("X_COORD_CD", "Y_COORD_CD"), crs = 2263)
  
  cdta <- st_read("./data/nycdta2020_24d/nycdta2020.shp")
  
  cdta <- st_transform(cdta, crs = st_crs(df_sf))
  
  df_joined <- st_join(df_sf, cdta, left = TRUE)
  
  df$CDTA <- df_joined$CDTA2020
  
  df$CDTA <- gsub("([A-Z]{2})(\d{2})", "\1 \2", df$CDTA)

• filters the poverty and education data for the relevant time period (2017-21) and geographic type (NTA2020), selects specific columns, and merges this poverty data with the shooting incident dataset.
  detailed code

  neighborhood_poverty <- read.csv("./data/neighborhood_poverty.csv")
  
  # Select specific columns from neighborhood_poverty
  neighborhood_poverty_selected <- neighborhood_poverty %>% 
    filter(TimePeriod == '2017-21') %>% 
    filter(GeoType == 'NTA2020' ) %>% 
    select(Number, Percent, Geography)
  # Merge the dataframes on 'NTA' and 'Geography', keeping all information from data_processing
  df_poverty <- left_join(df, neighborhood_poverty_selected, by = c("NTA" = "Geography"))
  
  # Rename columns after merging
  df_poverty <- df_poverty %>% 
    rename(Number_poverty = Number, Percent_poverty = Percent)
  neighborhood_education = read.csv("./data/graduated_high_school.csv")
  neighborhood_education_selected <- neighborhood_education %>% 
    filter(TimePeriod == '2017-21') %>% 
    filter(GeoType == 'NTA2020' ) %>% 
    select(Number, Percent, Geography)
  # Merge the dataframes on 'NTA' and 'Geography', keeping all information from data_processing
  df_education <- left_join(df_poverty, neighborhood_education_selected, by = c("NTA" = "Geography"))
  
  # Rename columns after merging
  df_education <- df_education %>% 
    rename(Number_education = Number, Percent_education = Percent)
  # filter between 2017 and 2023
  df_education$OCCUR_DATE <- as.Date(df_education$OCCUR_DATE, format = "%Y-%m-%d")
  
  # Filter the incidents that happened between 2017 and 2023
  df_education <- df_education %>% 
    filter(OCCUR_DATE >= as.Date("2017-01-01") & OCCUR_DATE <= as.Date("2023-12-31"))

• calculates shooting incident rates for both Neighborhood Tabulation Areas (NTAs) and boroughs by year.
  detailed code

  #add a column that shows the shooting incident rate of NTAs in that year
  df_education <- df_education %>%
    group_by(NTA, Year) %>%
    mutate(incident_rate_by_year_nta = (n() / Total_population)*100) %>%
    ungroup()
  
  #add a column that shows the shooting incident rate of boroughs in that year
  boro_population <- df_education %>%
    group_by(BORO, Year) %>%
    summarise(total_population_boro = sum(Total_population, na.rm = TRUE)) %>%
    ungroup()
  
  df_education <- df_education %>%
    left_join(boro_population, by = c("BORO", "Year")) %>%
    group_by(BORO, Year) %>%
    mutate(incident_rate_by_year_boro = (n() / total_population_boro)*100) %>%
    ungroup()
  
  #Rename a column
  df_education <- df_education %>%
    rename(Total_population_nta = Total_population)
  
  write.csv(df_education, "data_final.csv", row.names = FALSE)

The detailed code is here

Neighborhood Poverty

neighborhood_poverty <- read.csv("./data/neighborhood_poverty.csv")

# Select specific columns from neighborhood_poverty
neighborhood_poverty_selected <- neighborhood_poverty %>% 
  filter(TimePeriod == '2017-21') %>% 
  filter(GeoType == 'NTA2020' ) %>% 
  select(Number, Percent, Geography)

• Downloaded the full table from the NYC Environment and Health Data Portal
• Filtered to include only percent data from 2017-2021 for each NTA
• Describes estimated percentage of people whose annual income falls below 100% of the federal poverty level in each NTA

Graduated High School

neighborhood_education=read.csv("./data/graduated_high_school.csv")
neighborhood_education_selected <- neighborhood_education %>% 
  filter(TimePeriod == '2017-21') %>% 
  filter(GeoType == 'NTA2020' ) %>% 
  select(Number, Percent, Geography)

• Downloaded the full table from the NYC Environment and Health Data Portal
• Filtered the data to include only the percentage information for high school graduation between 2017 and 2021, focusing on Neighborhood Tabulation Areas (NTAs) as defined by the 2020 census.
• This dataset describes the estimated percentage of individuals aged 25 and older who completed high school or obtained a high school equivalency in each NTA

Exploratory data analysis

Maps

Statistical analysis

Method for Calculating Correlation Coefficients

To assess potential linear relationships between socioeconomic variables and shooting rates, correlation coefficients were calculated. Specifically, the relationship between poverty rates and shooting rates, as well as between the percentage of high school graduates and shooting rates, was explored across different neighborhoods in NYC.

Poverty

Calculate the correlation between the poverty percentage and the incident rate.

correlation <- cor(data_clean$incident_rate_by_year_nta, data_clean$Percent_poverty, use = "complete.obs")
print(paste("Correlation coefficient: ", correlation))

## [1] "Correlation coefficient:  0.508408410575774"

data_clean %>%
  plot_ly(x = ~Percent_poverty, y = ~incident_rate_by_year_nta, 
          color = ~NTA, colors = "viridis", 
          type = "scatter", mode = "markers",
          text = ~paste("Neighborhood: ", NTA, "<br>Borough: ", BORO, 
                        "<br>% Below Poverty Line: ", Percent_poverty, 
                        "<br>Incident Rate: ", incident_rate_by_year_nta)) %>%
  layout(title = "Percent Below the Poverty Line and Incident Rate in NYC",
         xaxis = list(title = 'Percentage of People Whose Income is Below the Poverty Line'),
         yaxis = list(title = 'Incident Rate'),
         legend = list(title = list(text = 'Neighborhood')))

## Warning: Ignoring 95 observations

# Scatter plot for Brooklyn
data_clean |> 
  filter(neighbourhood_group == "Brooklyn") |> 
  plot_ly(data = _, x = ~Percent_poverty, y = ~incident_rate_by_year_nta, 
          color = ~NTA,
          colors = "plasma", 
          type = "scatter",
          mode = "markers",
          text = ~paste("Neighborhood: ", NTA, "<br>Borough: ", neighbourhood_group, 
                        "<br>% Below Poverty Line: ", Percent_poverty, 
                        "<br>Incident Rate: ", incident_rate_by_year_nta)) |> 
    layout(title = "Percent Below the Poverty Line and Incident Rate in Brooklyn",
           xaxis = list(title = 'Percentage of People Whose Income is Below the Poverty Line'),
           yaxis = list(title = 'Incident Rate'),
           legend = list(title = list(text = 'Neighborhood')))

## Warning: Ignoring 13 observations

# Scatter plot for Staten Island
data_clean |> 
  filter(neighbourhood_group == "Staten Island") |> 
  plot_ly(data = _, x = ~Percent_poverty, y = ~incident_rate_by_year_nta, 
          color = ~NTA,
          colors = "inferno", 
          type = "scatter",
          mode = "markers",
          text = ~paste("Neighborhood: ", NTA, "<br>Borough: ", neighbourhood_group, 
                        "<br>% Below Poverty Line: ", Percent_poverty, 
                        "<br>Incident Rate: ", incident_rate_by_year_nta)) |> 
    layout(title = "Percent Below the Poverty Line and Incident Rate in Staten Island",
           xaxis = list(title = 'Percentage of People Whose Income is Below the Poverty Line'),
           yaxis = list(title = 'Incident Rate'),
           legend = list(title = list(text = 'Neighborhood')))

## Warning: Ignoring 5 observations

• Across all neighborhoods in NYC, there is a moderate positive linear relationship between the percentage of people below the poverty line (Percent_poverty) and the incident rate by neighborhood (incident_rate_by_year_nta).(r = 0.508).
• Associations differ by borough.
• Notably, Brooklyn (0.5686) and Staten Island (0.5576) show the highest correlations, suggesting that the relationship between poverty and incident rate is stronger in these boroughs.

Education

Calculate the correlation between the graduated in highschool percentage and the incident rate

# Calculate the correlation between the graduated in highschool percentage and the incident rate
correlation <- cor(data_clean$incident_rate_by_year_nta, data_clean$Percent_education, use = "complete.obs")
print(paste("Correlation coefficient: ", correlation))

## [1] "Correlation coefficient:  -0.274782643621427"

# Create a scatter plot to visualize the relationship
data_clean %>%
  plot_ly(x = ~Percent_education, y = ~incident_rate_by_year_nta, 
          color = ~NTA, colors = "viridis", 
          type = "scatter", mode = "markers",
          text = ~paste("Neighborhood: ", NTA, "<br>Borough: ", BORO, 
                        "<br>% graduated HS: ", Percent_education, 
                        "<br>Incident Rate: ", incident_rate_by_year_nta)) %>%
  layout(title = "Percent graduated high school and Incident Rate in NYC",
         xaxis = list(title = 'Percentage of People graduated in high school'),
         yaxis = list(title = 'Incident Rate'),
         legend = list(title = list(text = 'Neighborhood')))

## Warning: Ignoring 95 observations

# Scatter plot for The Bronx
data_clean |> 
  filter(neighbourhood_group == "Bronx") |>
  plot_ly(data = _, x = ~Percent_poverty, y = ~incident_rate_by_year_nta, 
          color = ~NTA,
          colors = "magma", 
          type = "scatter",
          mode = "markers",
          text = ~paste("Neighborhood: ", NTA, "<br>Borough: ", neighbourhood_group, 
                        "<br>% graduated HS: ", Percent_education, 
                        "<br>Incident Rate: ", incident_rate_by_year_nta)) |> 
   layout(title = "Percent graduated high school and Incident Rate in The Bronx",
           xaxis = list(title = 'Percentage of People graduated in high school'),
           yaxis = list(title = 'Incident Rate'),
           legend = list(title = list(text = 'Neighborhood')))

## Warning: Ignoring 25 observations

# Scatter plot for Staten Island
data_clean |> 
  filter(neighbourhood_group == "Staten Island") |>
  plot_ly(data = _, x = ~Percent_poverty, y = ~incident_rate_by_year_nta, 
          color = ~NTA,
          colors = "inferno", 
          type = "scatter",
          mode = "markers",
          text = ~paste("Neighborhood: ", NTA, "<br>Borough: ", neighbourhood_group, 
                        "<br>% graduated HS: ", Percent_education, 
                        "<br>Incident Rate: ", incident_rate_by_year_nta)) |> 
   layout(title = "Percent graduated high school and Incident Rate in Staten Island",
           xaxis = list(title = 'Percentage of People graduated in Staten Island'),
           yaxis = list(title = 'Incident Rate'),
           legend = list(title = list(text = 'Neighborhood')))

## Warning: Ignoring 5 observations

The correlation coefficient is -0.2748, indicating a negative relationship between the percentage of people who graduated from high school and the incident rate. This means that there is a tendency for higher education levels to be associated with lower incident rates * Across all neighborhoods in NYC,there is no clear linear trend between high school graduation rates and incident rates in Manhattan. * Associations differ by borough. * Notably, The Bronx (-0.3996) and Staten Island (-0.6452) show the highest correlations, suggesting that the relationship between percentage of people who graduated from high school and incident rate is stronger in these boroughs.

Discussion

Limitations